Abortion Access in the US¶

Hanqing Ye¶

Urban Economics 2023¶

This project looks at abortion access in the US over time. I obtain two variables from the Guttmacher Institute: Percentage of patients traveling out of state for abortion, and number of abortion providers in each state. For the final map, I make a GIF of the percentage of patients traveling out of state to obtain care from 2011 to 2020.

In [ ]:
import pandas as pd
import geopandas as gpd
import os
import xarray, rioxarray
import contextily
import seaborn as sns
from pysal.viz import mapclassify as mc
from legendgram import legendgram
import matplotlib.pyplot as plt
import palettable.matplotlib as palmpl
from splot.mapping import vba_choropleth
from shapely.geometry import Polygon
import warnings
warnings.filterwarnings('ignore')

os.chdir('/Users/hanqingye/Desktop/urban_map')

# Providers data
providers = pd.read_csv(os.getcwd()+'/data/no-of-abortion-providers.csv')
providers.head()
Out[ ]:
measure_name datum state_id state_name datum_date first_year last_year footnotes sources
0 No. of abortion providers 32 AL Alabama 1/1/78 1978.0 NaN A provider is a hospital, clinic or physician'... Abortion Incidence and Access to Services in t...
1 No. of abortion providers 23 AK Alaska 1/1/78 1978.0 NaN A provider is a hospital, clinic or physician'... Abortion Incidence and Access to Services in t...
2 No. of abortion providers 39 AZ Arizona 1/1/78 1978.0 NaN A provider is a hospital, clinic or physician'... Abortion Incidence and Access to Services in t...
3 No. of abortion providers 12 AR Arkansas 1/1/78 1978.0 NaN A provider is a hospital, clinic or physician'... Abortion Incidence and Access to Services in t...
4 No. of abortion providers 527 CA California 1/1/78 1978.0 NaN A provider is a hospital, clinic or physician'... Abortion Incidence and Access to Services in t...
In [ ]:
# Keep columns needed and change column names
col_names = ['datum', 'state_id', 'state_name', 'datum_date']
providers = providers[col_names]
new_col_names = ['num_providers', 'state_id', 'state_name', 'year']
providers.columns = new_col_names
providers['year'] = pd.to_datetime(providers['year']).dt.year
providers.dtypes
Out[ ]:
num_providers     int64
state_id         object
state_name       object
year              int32
dtype: object
In [ ]:
# Traveling out of state to obtain care data
travel_oos = pd.read_csv(os.getcwd()+'/data/of-residents-obtaining-abortions-who-traveled-out-of-state-for-care.csv')
travel_oos.head()
Out[ ]:
measure_name datum state_id state_name first_year last_year footnotes sources
0 % of residents obtaining abortions who travele... 17.0 AL Alabama 2011 NaN NaN 1.\tMaddow-Zimet I and Kost K, Even before Roe...
1 % of residents obtaining abortions who travele... 11.0 AK Alaska 2011 NaN NaN 1.\tMaddow-Zimet I and Kost K, Even before Roe...
2 % of residents obtaining abortions who travele... 2.0 AZ Arizona 2011 NaN NaN 1.\tMaddow-Zimet I and Kost K, Even before Roe...
3 % of residents obtaining abortions who travele... 27.0 AR Arkansas 2011 NaN NaN 1.\tMaddow-Zimet I and Kost K, Even before Roe...
4 % of residents obtaining abortions who travele... 0.0 CA California 2011 NaN NaN 1.\tMaddow-Zimet I and Kost K, Even before Roe...
In [ ]:
# Keep columns needed and change column names
col_names = ['datum', 'state_id', 'state_name', 'first_year']
travel_oos = travel_oos[col_names]
new_col_names = ['out_of_state_pct','state_id', 'state_name', 'year']
travel_oos.columns = new_col_names
travel_oos.dtypes
Out[ ]:
out_of_state_pct    float64
state_id             object
state_name           object
year                  int64
dtype: object
In [ ]:
# join two dataframes
abortion = pd.merge(providers, travel_oos, on=['state_id', 'state_name', 'year'], how='outer')
abortion.head()
Out[ ]:
num_providers state_id state_name year out_of_state_pct
0 32.0 AL Alabama 1978 NaN
1 23.0 AK Alaska 1978 NaN
2 39.0 AZ Arizona 1978 NaN
3 12.0 AR Arkansas 1978 NaN
4 527.0 CA California 1978 NaN

Geodata¶

In [ ]:
# US states geo data
gdf = gpd.read_file(os.getcwd()+'/data/cb_2018_us_state_500k')
gdf.head()
Out[ ]:
STATEFP STATENS AFFGEOID GEOID STUSPS NAME LSAD ALAND AWATER geometry
0 28 01779790 0400000US28 28 MS Mississippi 00 121533519481 3926919758 MULTIPOLYGON (((-88.50297 30.21523, -88.49176 ...
1 37 01027616 0400000US37 37 NC North Carolina 00 125923656064 13466071395 MULTIPOLYGON (((-75.72681 35.93584, -75.71827 ...
2 40 01102857 0400000US40 40 OK Oklahoma 00 177662925723 3374587997 POLYGON ((-103.00257 36.52659, -103.00219 36.6...
3 51 01779803 0400000US51 51 VA Virginia 00 102257717110 8528531774 MULTIPOLYGON (((-75.74241 37.80835, -75.74151 ...
4 54 01779805 0400000US54 54 WV West Virginia 00 62266474513 489028543 POLYGON ((-82.64320 38.16909, -82.64300 38.169...
In [ ]:
gdf.explore()
Out[ ]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [ ]:
# Join abortion data with geo data
df_final = pd.merge(gdf, abortion, left_on='STUSPS', right_on='state_id', how='outer')
df_final.head()
Out[ ]:
STATEFP STATENS AFFGEOID GEOID STUSPS NAME LSAD ALAND AWATER geometry num_providers state_id state_name year out_of_state_pct
0 28 01779790 0400000US28 28 MS Mississippi 00 1.215335e+11 3.926920e+09 MULTIPOLYGON (((-88.50297 30.21523, -88.49176 ... 10.0 MS Mississippi 1978.0 NaN
1 28 01779790 0400000US28 28 MS Mississippi 00 1.215335e+11 3.926920e+09 MULTIPOLYGON (((-88.50297 30.21523, -88.49176 ... 9.0 MS Mississippi 1979.0 NaN
2 28 01779790 0400000US28 28 MS Mississippi 00 1.215335e+11 3.926920e+09 MULTIPOLYGON (((-88.50297 30.21523, -88.49176 ... 9.0 MS Mississippi 1980.0 NaN
3 28 01779790 0400000US28 28 MS Mississippi 00 1.215335e+11 3.926920e+09 MULTIPOLYGON (((-88.50297 30.21523, -88.49176 ... 14.0 MS Mississippi 1981.0 NaN
4 28 01779790 0400000US28 28 MS Mississippi 00 1.215335e+11 3.926920e+09 MULTIPOLYGON (((-88.50297 30.21523, -88.49176 ... 13.0 MS Mississippi 1982.0 NaN

Map for one year: 2011¶

Percentage of patients traveling out of state for abortions¶

In [ ]:
# Distribution
df_final.loc[df_final['year'] == 2011, 'out_of_state_pct'].hist()
plt.show()
No description has been provided for this image
In [ ]:
# Variable of interest
var = 'out_of_state_pct'

# year
year = 2011
df_plot = df_final[df_final['year'] == year]

# Set the bounds for coloring
vmin, vmax = 0, 80 

# Color scheme
colormap = 'OrRd'

# Center mainland US
fig, ax = plt.subplots(1, figsize=(18, 14))
ax.set(xlim=(-2.5e6, 3e6), ylim=(-2.5e6, 1.25e6))  
ax.axis('off')

# The title
hfont = {'fontname':'sans-serif'}
ax.set_title(f'Abortion Access: Percentage of patients \n traveling out of state to obtain abortion in {year}', **hfont, fontdict={'fontsize': '30', 'fontweight' : '1'}, y = 0.9)

# Plot mainland counties with no PDMP
df_mainland = df_plot[(~df_plot['STUSPS'].isin(['AK','HI']))]
df_mainland_plot = df_mainland.to_crs({'init':'epsg:2163'})
df_mainland_plot.plot(var,linewidth=0.8,ax=ax,vmin=vmin,vmax=vmax,edgecolor='0.8',aspect=1,missing_kwds= dict(color = "lightgrey"),cmap = colormap)

# Iterate over geometries and add text annotations
for idx, row in df_mainland_plot.iterrows():
    # Check if the geometry is not None
    if row.geometry is not None:
        # Get the centroid of the geometry
        centroid = row.geometry.centroid
        # Place the state_id as a text label
        ax.text(centroid.x, centroid.y, row['state_id'], fontsize=8, ha='center', va='center', fontdict={'fontname': 'sans-serif', 'fontsize': '10'})

# Add color bars
cbax = fig.add_axes([0.83, 0.20, 0.03, 0.31]) 
legend = plt.cm.ScalarMappable(cmap=colormap,norm=plt.Normalize(vmin=vmin, vmax=vmax))
legend._A = []
fig.colorbar(legend, cax=cbax, ticks= [0, 20,40,60,80,+80])
cbax.set_yticklabels(['0','20','40','60','80','+80'])

# Move Alaska to the bottom left
left_al, bottom_al, width_al, height_al = [0.15, 0.0, 0.2, 0.5] 
ax2 = fig.add_axes([left_al, bottom_al, width_al, height_al])
ax2.axis('off')

# plot Alaska 
alaska = df_plot[(df_plot['state_id'] == 'AK')]
polygon = Polygon([(-170,50),(-170,72),(-140, 72),(-140,50)])
alaska.clip(polygon).plot(var,linewidth=0.8,ax=ax2,vmin=vmin,vmax=vmax, edgecolor='0.8',aspect=1,missing_kwds= dict(color = "lightgrey"),cmap = colormap)

# Add state label
for idx, row in alaska.iterrows():
    if row.geometry is not None:
        centroid = row.geometry.centroid
        ax2.text(centroid.x, centroid.y, row['state_id'], fontsize=8, ha='center', va='center', fontdict={'fontname': 'sans-serif', 'fontsize': '10'})

# Move Hawaii to the bottom left
left_hi, bottom_hi, width_hi, height_hi = [0.35, 0.1, 0.1, 0.3]
ax3 = fig.add_axes([left_hi, bottom_hi, width_hi, height_hi])  
ax3.axis('off')

# Plot Hawaii
hawaii = df_plot[(df_plot['state_id'] == 'HI')]
hipolygon = Polygon([(-160,0),(-160,90),(-120,90),(-120,0)])
hawaii.clip(hipolygon).plot(var,linewidth=0.8,ax=ax3,vmin=vmin,vmax=vmax, edgecolor='0.8',aspect=1,missing_kwds= dict(color = "lightgrey"),cmap = colormap)

# Add state label
for idx, row in hawaii.iterrows():
    if row.geometry is not None:
        centroid = row.geometry.centroid
        ax3.text(centroid.x, centroid.y, row['state_id'], fontsize=8, ha='center', va='center', fontdict={'fontname': 'sans-serif', 'fontsize': '10'})
No description has been provided for this image

Making a GIF¶

In [ ]:
# All years of this data
oos_years = travel_oos['year'].unique()

# Variable of interest
var = 'out_of_state_pct'

# Set the bounds for coloring
vmin, vmax = 0, 80 

# Color scheme
colormap = 'OrRd'

# Create the output directory
output_dir = os.path.join(os.getcwd(), 'output')
if not os.path.exists(output_dir):
    os.makedirs(output_dir)

# Loop through each year
for year in oos_years:
    df_plot = df_final[df_final['year'] == year]

    # Create a new figure for each year
    fig, ax = plt.subplots(1, figsize=(18, 14))
    ax.set(xlim=(-2.5e6, 3e6), ylim=(-2.5e6, 1.25e6))  
    ax.axis('off')

    # The title
    hfont = {'fontname':'sans-serif'}
    ax.set_title(f'Abortion Access: Percentage of patients \n traveling out of state to obtain abortion in {year}', 
                 **hfont, fontdict={'fontsize': '30', 'fontweight' : '1'}, y = 0.9)

    # Plot mainland counties with no PDMP
    df_mainland = df_plot[(~df_plot['STUSPS'].isin(['AK','HI']))]
    df_mainland_plot = df_mainland.to_crs({'init':'epsg:2163'})
    df_mainland_plot.plot(var,linewidth=0.8,ax=ax,vmin=vmin,vmax=vmax,edgecolor='0.8',aspect=1,missing_kwds= dict(color = "lightgrey"),cmap = colormap)

    # Iterate over geometries and add text annotations
    for idx, row in df_mainland_plot.iterrows():
    # Check if the geometry is not None
        if row.geometry is not None:
        # Get the centroid of the geometry
            centroid = row.geometry.centroid
        # Place the state_id as a text label
            ax.text(centroid.x, centroid.y, row['state_id'], fontsize=8, ha='center', va='center', fontdict={'fontname': 'sans-serif', 'fontsize': '10'})

    # Add color bars
    cbax = fig.add_axes([0.83, 0.20, 0.03, 0.31]) 
    legend = plt.cm.ScalarMappable(cmap=colormap,norm=plt.Normalize(vmin=vmin, vmax=vmax))
    legend._A = []
    fig.colorbar(legend, cax=cbax, ticks= [0, 20,40,60,80,+80])
    cbax.set_yticklabels(['0','20','40','60','80','+80'])

# Move Alaska to the bottom left
    left_al, bottom_al, width_al, height_al = [0.15, 0.0, 0.2, 0.5] 
    ax2 = fig.add_axes([left_al, bottom_al, width_al, height_al])
    ax2.axis('off')

# plot Alaska 
    alaska = df_plot[(df_plot['state_id'] == 'AK')]
    polygon = Polygon([(-170,50),(-170,72),(-140, 72),(-140,50)])
    alaska.clip(polygon).plot(var,linewidth=0.8,ax=ax2,vmin=vmin,vmax=vmax, edgecolor='0.8',aspect=1,missing_kwds= dict(color = "lightgrey"),cmap = colormap)

# Add state label
    for idx, row in alaska.iterrows():
        if row.geometry is not None:
            centroid = row.geometry.centroid
            ax2.text(centroid.x, centroid.y, row['state_id'], fontsize=8, ha='center', va='center', fontdict={'fontname': 'sans-serif', 'fontsize': '10'})

# Move Hawaii to the bottom left
    left_hi, bottom_hi, width_hi, height_hi = [0.35, 0.1, 0.1, 0.3]
    ax3 = fig.add_axes([left_hi, bottom_hi, width_hi, height_hi])  
    ax3.axis('off')

# Plot Hawaii
    hawaii = df_plot[(df_plot['state_id'] == 'HI')]
    hipolygon = Polygon([(-160,0),(-160,90),(-120,90),(-120,0)])
    hawaii.clip(hipolygon).plot(var,linewidth=0.8,ax=ax3,vmin=vmin,vmax=vmax, edgecolor='0.8',aspect=1,missing_kwds= dict(color = "lightgrey"),cmap = colormap)

# Add state label
    for idx, row in hawaii.iterrows():
        if row.geometry is not None:
            centroid = row.geometry.centroid
            ax3.text(centroid.x, centroid.y, row['state_id'], fontsize=8, ha='center', va='center', fontdict={'fontname': 'sans-serif', 'fontsize': '10'})


    # Save the figure with a unique name for each year
    plt.savefig(os.path.join(output_dir, f'abortion_access_{year}.png'), bbox_inches='tight')

    plt.close(fig)  # Close the figure to free memory
In [ ]:
# Make GIF
from PIL import Image
import os
from IPython.display import Image as DisplayImage, display

output_path = '/Users/hanqingye/Desktop/urban_map/output'  # Replace with your folder path
resized_path = '/Users/hanqingye/Desktop/urban_map/resized_output'  # Replace with your desired output folder path

# Desired size
width, height = 1400, 1100  # Replace with your desired size

# Create the output folder if it doesn't exist
if not os.path.exists(resized_path):
    os.makedirs(resized_path)

# Resize images
for filename in os.listdir(output_path):
    if filename.endswith('.png'):  # or '.jpg', etc
        img = Image.open(os.path.join(output_path, filename))
        img = img.resize((width, height), Image.Resampling.LANCZOS)  # Updated here
        img.save(os.path.join(resized_path, filename))


# Create a list to hold the images
images = []

# Load each file into a list
for filename in sorted(os.listdir(resized_path)):
    if filename.endswith('.png'):  # or '.jpg', etc
        file_path = os.path.join(resized_path, filename)
        img = Image.open(file_path) 
        images.append(img)

# Save the images as a GIF
gif = '/Users/hanqingye/Desktop/urban_map/output/travel_oos.gif' 

# Save the first image and append the rest as frames
images[0].save(gif, save_all=True, append_images=images[1:], optimize=False, duration=1000, loop=0)

# Display the GIF
display(DisplayImage(filename=gif))

<IPython.core.display.Image object>
In [ ]:
# Display GIF in HTML export
import base64
from IPython.display import HTML

# Path to your GIF file
gif_path = '/Users/hanqingye/Desktop/urban_map/output/travel_oos.gif'  # Update with your actual path

with open(gif_path, 'rb') as gif_file:
    gif_base64 = base64.b64encode(gif_file.read()).decode('utf-8')

# Embed the GIF in HTML
html_str = f'<img src="data:image/gif;base64,{gif_base64}" alt="GIF">'
display(HTML(html_str))
GIF

Number of providers¶

In [ ]:
df_final.loc[df_final['year'] == 2011, 'num_providers'].hist()
plt.show()
No description has been provided for this image
In [ ]:
# Variable of interest
var = 'num_providers'

# year
year = 2011
df_plot = df_final[df_final['year'] == year]

# Set the bounds for coloring
vmin, vmax = 0, 250 

# Color scheme
colormap = 'BuPu'

# Center mainland US
fig, ax = plt.subplots(1, figsize=(18, 14))
ax.set(xlim=(-2.5e6, 3e6), ylim=(-2.5e6, 1.25e6))  
ax.axis('off')

# The title
hfont = {'fontname':'sans-serif'}
ax.set_title(f'Abortion Access: \n Number of abortion providers in {year}', **hfont, fontdict={'fontsize': '30', 'fontweight' : '1'})

# Plot mainland counties with no PDMP
df_mainland = df_plot[(~df_plot['STUSPS'].isin(['AK','HI']))]
df_mainland_plot = df_mainland.to_crs({'init':'epsg:2163'})
df_mainland_plot.plot(var,linewidth=0.8,ax=ax,vmin=vmin,vmax=vmax,edgecolor='0.8',aspect=1,missing_kwds= dict(color = "lightgrey"),cmap = colormap)

# Iterate over geometries and add text annotations
for idx, row in df_mainland_plot.iterrows():
    # Check if the geometry is not None
    if row.geometry is not None:
        # Get the centroid of the geometry
        centroid = row.geometry.centroid
        # Place the state_id as a text label
        ax.text(centroid.x, centroid.y, row['state_id'], fontsize=8, ha='center', va='center')

cbax2 = fig.add_axes([0.92, 0.21, 0.03, 0.31]) 
legend_2 = plt.cm.ScalarMappable(cmap=colormap,norm=plt.Normalize(vmin=vmin, vmax=vmax))
legend_2._A = []
fig.colorbar(legend_2, cax=cbax2, ticks= [0, 50,100,150,200,+250])
cbax2.set_yticklabels(['0','50','100','150','200','+250'])
Out[ ]:
[Text(1, 0, '0'),
 Text(1, 50, '50'),
 Text(1, 100, '100'),
 Text(1, 150, '150'),
 Text(1, 200, '200'),
 Text(1, 250, '+250')]
No description has been provided for this image
In [ ]: